Kurzeinführung in die Datenanlyse mit R

Dies ist eine kurze Einführung in R für (Noch-)Nichtprogrammierer zum Einlesen, Analysieren und Darstellen exemplarischer Daten. Das Ziel ist nicht, Euch zu perfekten R ProgrammiererInnen zu machen, sondern durch Aufzeigen einiger praktischer und einfacher Funktionen die Verwendung von Programmiersprachen/Entwicklungsumgebungen schmackhaft zu machen und etwas “die Scheu” zu nehmen diese zu verwenden!

Neben der hier dargestellten Funktionalität könnt Ihr sehr einfach weitere Informationen zum Arbeiten mit R im Internet finden, z.B. durch Verwendung der Google Suchfunktion und den Ergebnissen auf der Seite stackoverflow.com - einfach mal ausprobieren.

Die hier verwendeten Beispieldaten sind Simulationen eines Energie-Bilanz-Schneemodels für den Standort Col de Porte (Frankreich) für die Saison 2005/2006 sowie Beobachtungen des Schneewasseräquivalents für den selben Standort und Zeitraum. Informationen zur Station und den Daten findet Ihr hier:

https://essd.copernicus.org/preprints/5/29/2012/essdd-5-29-2012.pdf

Die hier verwendeten Daten können hier heruntergeladen werden:

https://www.acwr.eu/Protected/R/Data/Data_ColDePorte.zip

1 - Arbeitsumgebungen

Das Arbeiten mit R erfordert zunächst das R-Basispaket, welches für die jeweilig genutzten Betriebssysteme heruntergeladen und installiert werden muss. Pakete für die unterschiedlichen Betriebssysteme gibt es hier: https://cran.r-project.org/

Auf die dabei installierte Funktionalität (die R-Skriptsprache) kann man über unterschiedliche Umgebungen zugreifen:

R in der Konsole

Alle Betriebssysteme verfügen über eine “Konsole”, in dieser kann man verscheidene Kommandos eingeben um z.B. Verzeichnisse zu wechseln.

Befehl R: direktes Ausführen von R-Befehlen in der Konsole

Befehl Rscript myscript.R: Ausführen bestehender Skripte, hier “myscript.R”

R im Jupyter Notebook

Das Jupyter Notebook stellt eine übersichtliche Oberfläche zum Arbeiten in R dar, ist aber was die Datenverabreitung und den Funktionsumfang angeht in mancher Hinsicht etwas limitiert.

RStudio

RStudio ist die wohl am meisten verwendete Umgebung: https://www.rstudio.com/products/rstudio/download/

Hier kann man Code schreiben, Pakete installieren, Variablen bequem einsehen und Plots betrachten. Wir arbeiten in R-Studio, macht es doch gleich mal auf und öffnet ein neues Skript “Introduction.R” das Ihr im Ordner “R-Workspace” auf Eurem Rechner speichert.

Die Daten der Station “Col de Porte” (Frankreich) ladet Ihr auch gleich herunter und legt diese im gleichen Ordner ab!

2 - Grundlegendes

Programmablauf

Ein R-Script wird Zeile für Zeile von oben nach unten ausgeführt. So können Variablen (z.B. Name) mit Werten (z.B. Alex) gefüllt werden und für numerische Variablen Berechnungen durchgeführt werden.

Kommentare

In jeder Programmiersprache ist es wichtig den Code (die Befehlsabfolge) von Kommentaren (Erläuterungen die ein Lesen des Codes vereinfachen und den Inhalt für Andere nachvollziehbar machen) zu trennen. Kommentare können im R-Code durch # eingeleitet werden - diese Teile werden dann nicht ausgeführt, sondern sind nur für den Entwickler von Wert. Uns ist das Kommentieren sehr wichtig, da es den LV-Leitern ermöglicht Euren Code zu lesen und zu verstehen - bitte kommentiert doch deshalb Euren Code zur besseren Nachvollziehbarkeit - das hilft auch dabei vor langer Zeit geschriebene Skripte wieder verstehen/verwenden zu können.

Unten sehr Ihr ein Beispiel für einen Kommentar:

# Hello ausgeben
print('Hello')
## [1] "Hello"

Ausführen von Code

Code kann in R zeilenweise, für Bereiche oder für das gesamte Skript ausgeführt werden - hierfür gibt es fest definierte Shortcuts:

Ausführen der aktutellen Zeile/eines markierten Bereiches: CTRL & Enter

Ausführen des gesamten Skriptes: Shift & CTRL & s

Variablen

In der Programmierung arbeitet man mit “Variablen” die durch “Operationen” miteinander verarbeitet werden. So kann man einer frei zu definierenden Variablen einen Wert zuweisen und diese Variable dann mit einer anderen verrechnen.

Die Wertzuweisung erfolg in R durch <- oder =

# Erste Berechnungen
apples = 11
oranges = 4
fruits = apples + oranges

Schreiben wir den Namen der Variablen in eine Zeile und führen diese Zeile aus wird der Wert der Variable ausgegeben:

# Wert ausgeben
fruits
## [1] 15

Wir können uns den Typ einer Variablen durch class() anzeigen lassen.

In diesem Fall ist es eine Variable vom Typ numeric, also eine Zahl. Neben Zahlen gibt es zum Beispiel noch character Variablen, diese enthalten eine Zeichenfolge. “15” könnte also sowohl eine numerische Variable sein, aber auch eine Zeichenfolge - mit Zeichenfolgen kann man aber natürlich nicht rechnen, es empfiehlt sich bei Problemen im workflow immer mal den Typ einer Variablen (v.a. bei Zahlen) zu checken, das erspart oft viel Kopfzerbrechen:

# Klasse ausgeben
class(fruits)
## [1] "numeric"

Der Befehl str() gibt uns Information über die Struktur unserer Daten.

In diesem Fall ist das nur die Information, das es sich um eine Zahl mit dem Wert 15 handelt. Wir werden später mit DataFrames arbeiten, hier können unterschiedliche Spalten unterschiedliche Typen (z.B. Zahl, Zeichenfolge) haben.

# Struktur ausgeben
str(fruits)
##  num 15

Einbinden von Bibliotheken (Libraries)

Während viele Befehle zur Standardfunktion von R gehören, erfordern einige Operationen Zusatzfunktionalität, welche in Form von Bibliotheken eingebunden werden kann.

Entsprechende Bibliotheken können in RStudio im Fenster rechts unten durch Packages-->Install installiert werden, oder direkt als R-Code durch den Befehl install.packages("Paketname").

Wir installieren als Beispiel die Bibliothek this.path, die uns später den Umgang mit Arbeitsverzeichnissen erleichtert

Um das Paket im Code verwenden zu können, muss der Befehl library() in Zusammenhang mit dem Namen der einzubindenden Bibliothek eingeben werden. Wir binden Testweise die Library this.path ein:

library("this.path")

In R arbeitet man oft mit vordefinierten Funktionen die in Bibliotheken definiert sind. Will ich Hilfe zu einer solchen Funktion, hilft Aufruf von help(Funktionsname). Wir testen das anhand der Mittelwertsfunktion mean():

help(mean)

Das Arbeitsverzeichnis

Habe ich irgendwo ein R-Skript angelegt, so ist es oft praktisch das Arbreitsverzeichnis auf den Ordner mit diesem Skript zu setzen. Da Daten in R standardmäßig (soweit nicht mit absolutem Pfad definiert) im Arbeitsverzeichnis gesucht werden, kann ich auf Daten die direkt im Arbeitsverzeichnis leigen ohne Angabe eines Pfads nur über den Dateinamen zugreifen. Liegen Daten in einem Unterverzeichnis des Arbeitsverzeichnisses, dann muss ich diesen Ordner relativ zum Arbeitsverzeichnis schon noch angeben.

Das Arbeitsverzeichnis kann durch den Befehl setwd("/MeinOrdner") im R-Code gesetzt werden. Alternativ kann in RStudio der aktuell angezeigte Ordner im Fenster rechts unten als Arbeitsverzeichnis gesetzt werden Files-->More-->Set As Working Directory.

Durch den Befehl getwd() kann ich mir das aktuelle Arbeitsverzeichnis anzeigen lassen.

getwd()
## [1] "/Users/t/pCloud Drive/Current/WS_2023/EU4 Gebirgsforschung/Introduction-R"

Die eben installierte Bibliothek this.path ist in diesem Zusammenhang besonders praktisch - sie ermöglicht es den Pfad in einem R-Skript automatisch direkt auf den Pfad des Skriptes zu setzen. So muss ich mir über die Lage meines R-Skriptes/Ordners auf meinem Rechner keine Gedanken machen, solange alle benötigten Daten im gleichen Ordner liegen.

# Bibliothek laden
library("this.path")

# Pfad des Skriptes auslesen und als Arbeitsverzeichnis setzen
script_path = this.path()
setwd(dirname(script_path))

3 - Einlesen von Daten

Das Einlesen tabellarischer Daten in R ist sehr einfach. Durch die Funktion read.table() können tabellarische Daten in unterschiedlichen Formaten eingelesen werden. Man gibt hierfür den Namen in ““, das Trennungszeichen (z.B. sep=",") sowie die Information darüber, ob ein Header (eine vorangehende Nichtdatenzeile mit Spaltennamen) existiert an (z.B. header=TRUE). Weiter kann ich angeben, welches Dezimaltrennzeichen in meinen Daten verwendet wird (z.B. dec = "."). Mit dem Argument skip= könnte man noch angeben, dass eine bestimmte Anzahl an Zeilen übersprungen werden soll, z.B. wenn in den ersten Zeilen unnötige Zusatzinformation enthalten ist.

Wir testen das mal mit einer Datei die Meteorologische Modellinputs sowie eine Vielzahl von Modellautputs enthält. Die Datei sollte bereits im Arbeitsverzeichnis liegen und heisst ModelOutput.csv.

Tipp:

Schaut Euch tabellarischen Daten immer zunächst im Texteditor Eurer Wahl (Windows: Notepad, Mac: TextEdit, oder einfach in RStudio) an. So seht Ihr wie die Daten aussehen (Spalten, Header, Trennzeichen, etc.)!

Optionen für Trennzeichen sep="":

# Daten einlesen
input_file = "ModelOutput.csv"
input_data = read.table(input_file, sep=",", dec = ".", header=TRUE)

Datenstruktur

Die eingelesenen Daten befinden sich nun im Speicher des Rechners als Werte der Variaben input_data und wir können uns diese Daten einmal ansehen. Der Befehl names(input_data) zeigt uns alle Spaltennamen an!

Tipp:

Unter Verwendung der Funktion setNames(NameDataFrame,c("Spaltenname1","Spaltename2",...)) kann man die Spalten jederzeit umbenennen (machen wir später, hier aber gerade nicht nötig).

# Spaltennamen ausgeben
names(input_data)
##  [1] "X"           "Year"        "Month"       "Day"         "Hour"       
##  [6] "ShortIn"     "LongIn"      "SolPrecip"   "LiqPrecip"   "AirTemp"    
## [11] "RelHum"      "WindSpeed"   "SurfPress"   "SWE"         "SnowAlbedo" 
## [16] "ShortOut"    "LongOut"     "SensFlux"    "LatFlux"     "ColdContent"
## [21] "Melt"        "Runoff"      "SnowDepth"   "LiqWatCont"  "IceCont"    
## [26] "SnowTemp"    "SnowDensity"

Frage ich die Klasse der Variablen input_data ab, so wird diese als data.frame angezeigt. DataFrames sind sehr ähnlich einer Excel Tabelle mit unterschiedlichen Zeilen und Spalten. Neben den DataFrames gibt es viele andere Formate (z.B. das xts Format) die alle Ihre Vorteile und evtl. leicht andere Befehlssyntax haben - aber bleiben wir bei unserem DataFrame…

# Klasse ausgeben
class(input_data)
## [1] "data.frame"

Den gesamten Inhalt der Daten kann ich mir leicht ansehen, ich muss dazu nur alle Werte der Variablen ausgeben. Das mache ich indem ich die Variable ohne weiteren Befehl ausführe.

Wir könnten diese Zeile nach dem Anschauen unserer Daten im Code auch mit # auskommentieren, damit der Output unseres Skriptes nicht unnötig unübersichtlich wird :-)

# Alle Daten ausgeben
#input_data

Will ich mir nur die ersten 10 Zeilen n=10 ansehen, kann ich die Funktion head() einsetzen! Wir sehen hier einen Fehlwert (“Not Available” = NA) in der Zeitreihe.

# Header (erste 10 Zeilen) ausgeben
head(input_data, n=10)
##     X Year Month Day Hour ShortIn LongIn SolPrecip LiqPrecip AirTemp RelHum
## 1   1 2005    10   1    0     0.0  283.1         0         0   277.8   78.2
## 2   2 2005    10   1    1     0.0  284.7         0         0   278.0   73.1
## 3   3 2005    10   1    2     0.0  285.8         0         0   277.7   76.1
## 4   4 2005    10   1    3     0.0  288.1         0         0   278.3   72.0
## 5   5 2005    10   1    4     0.0  293.9         0         0   277.7   73.0
## 6   6 2005    10   1    5     0.0  335.0         0         0   279.4   69.0
## 7   7 2005    10   1    6     1.4  347.2         0         0   280.7   57.8
## 8   8 2005    10   1    7    23.6  354.7         0         0   281.6   56.8
## 9   9 2005    10   1    8    96.4  346.0         0         0   282.9   59.8
## 10 10 2005    10   1    9   131.3  358.3         0         0   285.2   44.6
##    WindSpeed SurfPress SWE SnowAlbedo ShortOut LongOut SensFlux LatFlux
## 1        0.6     87480  NA         NA       NA      NA       NA      NA
## 2        0.0     87430  NA         NA       NA      NA       NA      NA
## 3        1.0     87390  NA         NA       NA      NA       NA      NA
## 4        0.5     87380  NA         NA       NA      NA       NA      NA
## 5        0.2     87360  NA         NA       NA      NA       NA      NA
## 6        1.2     87370  NA         NA       NA      NA       NA      NA
## 7        1.1     87330  NA         NA       NA      NA       NA      NA
## 8        1.1     87310  NA         NA       NA      NA       NA      NA
## 9        0.6     87320  NA         NA       NA      NA       NA      NA
## 10       0.6     87310  NA         NA       NA      NA       NA      NA
##    ColdContent Melt Runoff SnowDepth LiqWatCont IceCont SnowTemp SnowDensity
## 1           NA   NA     NA        NA         NA      NA       NA          NA
## 2           NA   NA     NA        NA         NA      NA       NA          NA
## 3           NA   NA     NA        NA         NA      NA       NA          NA
## 4           NA   NA     NA        NA         NA      NA       NA          NA
## 5           NA   NA     NA        NA         NA      NA       NA          NA
## 6           NA   NA     NA        NA         NA      NA       NA          NA
## 7           NA   NA     NA        NA         NA      NA       NA          NA
## 8           NA   NA     NA        NA         NA      NA       NA          NA
## 9           NA   NA     NA        NA         NA      NA       NA          NA
## 10          NA   NA     NA        NA         NA      NA       NA          NA

Ich kann die einzelnen Spalten der Variable ansprechen, indem ich das $ Zeichen verwende und den Spaltennamen angebe! Wir testen das mit einem einfachen Ausdrucken der Werte für die Stunde, wieder nur für die ersten 10 Zeilen!

# Header einer Spalte ausgeben
head(input_data$Hour, n=10)
##  [1] 0 1 2 3 4 5 6 7 8 9

Datenstatistik

Oft interessiert uns die Statistik der Daten, z.B. der Mittelwert oder die Standardabweichung - diese Information kann man sich über den Befehl summary ausgeben lassen! Wir sehen wieder die Information über enthaltene Fehlwerte NA in der Datenreihe.

# Zusammenfassung ausgeben
summary(input_data)
##        X             Year          Month             Day             Hour      
##  Min.   :   1   Min.   :2005   Min.   : 1.000   Min.   : 1.00   Min.   : 0.00  
##  1st Qu.:1639   1st Qu.:2005   1st Qu.: 3.000   1st Qu.: 8.00   1st Qu.: 5.75  
##  Median :3276   Median :2006   Median : 5.000   Median :16.00   Median :11.50  
##  Mean   :3276   Mean   :2006   Mean   : 6.033   Mean   :15.68   Mean   :11.50  
##  3rd Qu.:4914   3rd Qu.:2006   3rd Qu.:10.000   3rd Qu.:23.00   3rd Qu.:17.25  
##  Max.   :6552   Max.   :2006   Max.   :12.000   Max.   :31.00   Max.   :23.00  
##                                                                                
##     ShortIn          LongIn        SolPrecip        LiqPrecip       
##  Min.   :  0.0   Min.   :177.8   Min.   :0.0000   Min.   : 0.00000  
##  1st Qu.:  0.0   1st Qu.:268.1   1st Qu.:0.0000   1st Qu.: 0.00000  
##  Median :  0.0   Median :297.4   Median :0.0000   Median : 0.00000  
##  Mean   :102.3   Mean   :292.0   Mean   :0.0772   Mean   : 0.05946  
##  3rd Qu.:112.7   3rd Qu.:320.3   3rd Qu.:0.0000   3rd Qu.: 0.00000  
##  Max.   :995.1   Max.   :403.2   Max.   :9.1080   Max.   :10.26000  
##                                                                     
##     AirTemp          RelHum         WindSpeed        SurfPress    
##  Min.   :258.3   Min.   : 19.00   Min.   :0.0000   Min.   :84490  
##  1st Qu.:270.1   1st Qu.: 64.80   1st Qu.:0.1000   1st Qu.:86460  
##  Median :275.7   Median : 81.30   Median :0.6000   Median :87000  
##  Mean   :276.3   Mean   : 76.86   Mean   :0.9216   Mean   :86856  
##  3rd Qu.:281.7   3rd Qu.: 91.10   3rd Qu.:1.4000   3rd Qu.:87370  
##  Max.   :297.0   Max.   :102.20   Max.   :8.5000   Max.   :88410  
##                                                                   
##       SWE           SnowAlbedo        ShortOut          LongOut      
##  Min.   :  0.02   Min.   :0.5000   Min.   :-741.89   Min.   :-315.6  
##  1st Qu.:148.58   1st Qu.:0.7385   1st Qu.: -60.54   1st Qu.:-314.1  
##  Median :287.95   Median :0.8586   Median :   0.00   Median :-297.9  
##  Mean   :272.46   Mean   :0.8240   Mean   : -56.02   Mean   :-295.5  
##  3rd Qu.:392.83   3rd Qu.:0.9281   3rd Qu.:   0.00   3rd Qu.:-281.2  
##  Max.   :499.71   Max.   :0.9500   Max.   :   0.00   Max.   :-237.8  
##  NA's   :2615     NA's   :2615     NA's   :2615      NA's   :2615    
##     SensFlux          LatFlux          ColdContent            Melt       
##  Min.   :-67.595   Min.   :-57.9717   Min.   :-40.7333   Min.   :0.0000  
##  1st Qu.:  3.812   1st Qu.: -4.1091   1st Qu.:-11.0415   1st Qu.:0.0000  
##  Median : 13.787   Median :  0.0081   Median : -4.6945   Median :0.0000  
##  Mean   : 16.071   Mean   :  0.3702   Mean   : -7.1116   Mean   :0.1668  
##  3rd Qu.: 26.453   3rd Qu.:  5.2245   3rd Qu.: -0.3866   3rd Qu.:0.0000  
##  Max.   :119.472   Max.   : 61.5705   Max.   :  0.0000   Max.   :4.6933  
##  NA's   :2615      NA's   :2615       NA's   :2615       NA's   :2615    
##      Runoff         SnowDepth        LiqWatCont        IceCont     
##  Min.   :0.0000   Min.   :0.0002   Min.   : 0.000   Min.   :  0.0  
##  1st Qu.:0.0000   1st Qu.:0.6069   1st Qu.: 0.000   1st Qu.:148.5  
##  Median :0.0000   Median :0.9808   Median : 0.000   Median :287.9  
##  Mean   :0.1591   Mean   :0.9455   Mean   : 8.522   Mean   :263.9  
##  3rd Qu.:0.0000   3rd Qu.:1.2126   3rd Qu.:13.344   3rd Qu.:365.0  
##  Max.   :5.0213   Max.   :1.9417   Max.   :49.973   Max.   :488.0  
##  NA's   :2615     NA's   :2615     NA's   :2615     NA's   :2615   
##     SnowTemp      SnowDensity    
##  Min.   :254.5   Min.   : 10.74  
##  1st Qu.:265.4   1st Qu.:241.45  
##  Median :269.2   Median :270.24  
##  Mean   :268.6   Mean   :283.67  
##  3rd Qu.:272.8   3rd Qu.:322.92  
##  Max.   :273.1   Max.   :473.33  
##  NA's   :2615    NA's   :2615

4 - Linien und Balkendiagramme

Einfache Diagramme (= Plots) zu erstellen ist in R relativ einfach. Es gibt viele Plotting Bibliotheken (z.B. ggplot2), wir wollen mit der Standard Plottingfunktion plot() arbeiten, da diese für Einsteiger leichter zu verstehen ist.

Für Fortgeschrittene ist das Plotten mit ggplot2 eine tolle Alternative, da die Plots hier oft ohne viel Aufwand schöner werden, im Internet gibt es hierzu tolle Tutorials, einfach mal googeln :-)

Wir erstellen nun unseren ersten Plot und überlassen so gut wie alles der plot() Funktion, d.h. wir geben ausser den X-Werten und Y-Werten keine Argumente mit an die Funktion. Testweise plotten wir die Lufttemperatur AirTemp gegen den Index X:

# Lufttemperatur plotten
plot(input_data$X, input_data$AirTemp)

Das funktioniert schonmal, allerdings hat dieser Plot noch einige Mängel die wir nun anpassen wollen.

Um unsere Zeitreihe auf der X-Achse gut darstellen zu können, generieren wir nun eine Zeitinformation in unserem DataFrame (wir nennen die neue Spalte DateTime) aus den Spalten für Jahr, Monat, Tag und Stunde. Um diese Zeitinformation zusammenzustellen verwenden wir die R-Funktion paste() - sie ermöglicht das kombinieren von Text und Variablenwerten.

Wichtig:

Variablen werden hier durch Komma getrennt angegeben, dabei ist reiner Text in Anführungszeichen zu setzen, ganz am Ende wird nach einem abschliessenden Komma mit sep='' definiert wie die Einzeilteile zu trennen sind (hier ohne Inhalt weil wir keine automatsiche Trennung, z.B. durch Leerzeichen wollen, diese fügen wir lieber selbst ein).

# Spalte DateTime erzeugen
input_data$DateTime = paste(input_data$Year,'-',input_data$Month,'-',input_data$Day,' ',input_data$Hour,sep='')

Wir sehen uns nun mit head() und str() einmal an was in der Spalte DateTime abgelegt wurde:

# Header ausgeben
head(input_data$DateTime, n=10)
##  [1] "2005-10-1 0" "2005-10-1 1" "2005-10-1 2" "2005-10-1 3" "2005-10-1 4"
##  [6] "2005-10-1 5" "2005-10-1 6" "2005-10-1 7" "2005-10-1 8" "2005-10-1 9"
# Struktur ausgeben
str(input_data$DateTime)
##  chr [1:6552] "2005-10-1 0" "2005-10-1 1" "2005-10-1 2" "2005-10-1 3" ...

Noch ist unsere Zeitinformation eine Zeichenfolge, wir transformieren diese nun in eine POSIX Zeitinformation, dabei geben wir das Format mit, in dem unsere Zeichenfolge vorliegt, mit der Funktion head() geben wir uns gleich die ersten Zeilen aus:

# In Zeitinformation transformieren und Header ausgeben
input_data$DateTime = as.POSIXct(input_data$DateTime,format="%Y-%m-%d %H")
head(input_data$DateTime, n=10)
##  [1] "2005-10-01 00:00:00 CEST" "2005-10-01 01:00:00 CEST"
##  [3] "2005-10-01 02:00:00 CEST" "2005-10-01 03:00:00 CEST"
##  [5] "2005-10-01 04:00:00 CEST" "2005-10-01 05:00:00 CEST"
##  [7] "2005-10-01 06:00:00 CEST" "2005-10-01 07:00:00 CEST"
##  [9] "2005-10-01 08:00:00 CEST" "2005-10-01 09:00:00 CEST"

Nun plotten wir unsere Zeitreihe nocheinmal, diesmal die Lufttemperatur gegen unsere neue Zeitinformation:

# Lufttemperatur plotten
plot(input_data$DateTime, input_data$AirTemp)

Wir sehen, dass die plot-Funktion mit der Zeitinformation gut umgehen kann, ohne weitere Eingriffe werden in der X-Achse die Monate dargestellt.

Um weitere Mängel zu korrigieren plotten wir die Daten erneut - diesmal geben der plot() Funktion noch ein paar Informationen mehr mit, z.B.

  1. dass wir nur Linien, keine Punkte darstellen wollen –> type="l",

  2. dass wir gerne eine rote Linie hätten –> col="red" und

  3. dass wir die Beschriftung der Achsen gerne selbst übernehmen –> ann=FALSE.

# Lufttemperatur plotten
plot(input_data$DateTime, input_data$AirTemp, type="l", col="red", ann=FALSE)

Weitere Optionen für die Darsttellung der Daten type="":

Weiter Farboptionen col="":

Einen schönen Überblick über die in R standardmäßig (es gibt viele Erweiterungspakete) verfügbaren Farben findet Ihr unter: http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf

Nun wollen wir unseren Plot wie folgt erweitern:

  1. als zweite Zeitreihe soll die simulierte Schneetemperatur input_data$SnowTemp hinzugefügt werden

  2. ein *Titel, Beschriftungen für die X- und Y-Achse sowie eine Legende soll hinzugefügt werden

  3. wir wollen ein Grid hinter die Zeitreihen legen

Tipp:

Unterschiedliche Linientypen lty="":

Diese können auch über Nummern angesprochen werden, wollt Ihr mehr über linetypes wissen besucht: http://www.sthda.com/english/wiki/line-types-in-r-lty

# Unser ursprünglicher Plot von oben, erweitert um die Angabe des Linientyps
plot(input_data$DateTime, input_data$AirTemp, type="l", col="red", ann=FALSE, lty='solid')

# Zweite Zeitreihe SnowTemp hinzufügen
lines(input_data$DateTime, input_data$SnowTemp, type="l", col="blue", lty='twodash')

# Diagrammtitel hinzufügen
title(main="Lufttemperatur und simulierte Schneetemperatur \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Monat")
title(ylab="Temperatur (K)")

# Grid hinzufügen
grid(lty = 2, col = "gray", lwd = 1)

# Legende hinzufügen, leichten Versatz angeben
legend("topleft", inset=.05, c("Lufttemperatur","Schneetemperatur"), col=c("red","blue"), lty=c('solid','twodash'))

Wir wollen nun noch den Niederschlag als Balkendiagramm in blau plotten. Da wir am Gesamtniederschlag interessiert sind und in unseren Daten fester SolPrecip und flüssiger Niederschlag LiqPrecip getrennt vorliegen, addieren wir diese zunächst in einer neuen Spalte Precip. Dann plotten wir den Niederschlag unter Verwednung des Darstellungstypes h (type="h").

# Flüssigen und festen Niederschlag in neuer Spalte zu Gesamtniederschlag addieren
input_data$Precip = input_data$LiqPrecip + input_data$SolPrecip

# Niederschlag plotten
plot(input_data$DateTime, input_data$Precip, type="h", col="blue", ann=FALSE)

# Diagrammtitel hinzufügen
title(main="Niederschlag \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Monat")
title(ylab="Stündlicher Niederschlag (mm)")

Als nächstes wollen wir einen Dichteplot erzeugen, der uns zeigt, wie der Niederschlag über den Wertebereich verteilt ist. Hierzu verwenden wir erneut die plot() Funktion, berechnen uns aber mit der Funktion density() zunächst die Dichteverteilung density_Precip, die wir einfach einmal ausgeben um uns die Berechnungen einmal anzusehen:

# Dichtefunktion berechnen und ausgeben
density_Precip = density(input_data$Precip)
density_Precip
## 
## Call:
##  density.default(x = input_data$Precip)
## 
## Data: input_data$Precip (6552 obs.); Bandwidth 'bw' = 0.09585
## 
##        x                 y           
##  Min.   :-0.2875   Min.   :0.000000  
##  1st Qu.: 2.4212   1st Qu.:0.000341  
##  Median : 5.1300   Median :0.001970  
##  Mean   : 5.1300   Mean   :0.092235  
##  3rd Qu.: 7.8388   3rd Qu.:0.009271  
##  Max.   :10.5475   Max.   :3.677512

Nun wollen wir uns die berechnete Dichtefunktion density_Precip im Plot durch eine blaue Linie (col="blue") darstellen, wieder wollen wir Titel und Achsenbeschriftung selbst definieren (ann=FALSE).

# Plot erstellen
plot(density_Precip, col="blue", ann=FALSE)

# Diagrammtitel hinzufügen
title(main="Dichtefunktion des Niederschlags \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Niederschlag (mm)")
title(ylab="Dichte")

Wir fügen nun auch den festen Niederschlag hinzu, diesmal in einer exotischeren Farbe, darkturquoise

# Ursprünglicher Plot für den festen Niederschlag
plot(density_Precip, col="blue", ann=FALSE)

# Festen Niederschlag hinzufügen, vorher Dichteverteilung berechnen
density_SolPrecip <- density(input_data$SolPrecip)
lines(density_SolPrecip, col="darkturquoise")

# Diagrammtitel hinzufügen
title(main="Dichtefunktion des Niederschlags \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Niederschlag (mm)")
title(ylab="Dichte")

Da ist die zweite Linie wohl etwas über das Ziel hinaus geschossen. Das liegt daran, dass unser erster plot() Aufruf die Y-Werte auf den kleineren Wertebereich optimiert hat. Also versuchen wir es nochmal, diesmal begrenzen wir den Wertebereich der X- und Y-Achse durch die Argumete xlim und ylim. Die Angaben erfolgen hier in Form eines Vektors c() mit den Minimal- und Maximalwerten. Grid und Legende fügen wir auch gleich hinzu (beachtet, dass unser Grid zuerst hinzugefügt wird, sonst liegt es über der Legende):

# Ursprünglicher Plot für den festen Niederschlag
plot(density_Precip, xlim=c(0,12), ylim=c(0,6), col="blue", ann=FALSE)

# Festen Niederschlag hinzufügen, vorher Dichteverteilung berechnen
density_SolPrecip <- density(input_data$SolPrecip)
lines(density_SolPrecip, col="darkturquoise")

# Diagrammtitel hinzufügen
title(main="Dichtefunktion des Niederschlags \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Niederschlag (mm)")
title(ylab="Dichte")

# Grid hinzufügen
grid(lty = 'dotted', col = "gray", lwd = 1)

# Legende hinzufügen
legend("topright", inset=.05, c("Niederschlag","Flüssiger Niederschlag"), col=c("blue","darkturquoise"), lty=c("solid","solid"))

5 - Subsets aus den Daten auswählen

Manchmal will man nur mit einem **Zeitbereich in *einer Datenreihe** arbeiten. Dann kann man auch Zeiträume extrahieren und am besten einer anderen Variablen zuordnen - so verlieren wir nicht unsere ursprünglichen Daten.

Das Auswählen von Daten aus einer Datenreiche kann über unterschiedliche Wege erfolgen:

  1. Über den Index (z.B. Zeile 5-20)

  2. Über eine Abfrage der Inhalte (z.B. Datum X - Y)

  3. Über den Namen einer Spalte

Wir testen die Ansätze an einigen Beispielen unter Verwendung unseres eingelesenen Datensatzes input_data.

Zunächst wollen wir alle Spalten des DataFrames, aber nur die Zeilen 25 - 48 auswählen und in einen neuen DataFrame index_subset schreiben, wir schauen diesen gleich mit head() an:

# Bereich über den Index ausschneiden und Header ausgeben
index_subset = input_data[25:48,]
head(index_subset, n=50)
##     X Year Month Day Hour ShortIn LongIn SolPrecip LiqPrecip AirTemp RelHum
## 25 25 2005    10   2    0     0.0  339.7     0.000    0.4032   277.5   97.1
## 26 26 2005    10   2    1     0.0  322.8     0.000    0.0000   277.1   97.0
## 27 27 2005    10   2    2     0.0  331.4     0.000    0.0000   277.1   97.0
## 28 28 2005    10   2    3     0.0  327.2     0.000    0.0000   276.9   94.9
## 29 29 2005    10   2    4     0.0  326.7     0.000    0.0000   276.7   93.9
## 30 30 2005    10   2    5     0.0  328.6     0.000    0.7020   276.1   88.8
## 31 31 2005    10   2    6     0.3  325.0     0.000    0.9936   276.1   91.8
## 32 32 2005    10   2    7     3.9  328.2     0.000    3.1392   274.7   94.8
## 33 33 2005    10   2    8    14.4  317.7     0.000    2.4084   273.8   96.8
## 34 34 2005    10   2    9    35.3  324.4     0.000    1.4688   274.9   97.9
## 35 35 2005    10   2   10    53.1  327.2     0.000    3.4128   274.2   96.9
## 36 36 2005    10   2   11    43.3  314.4     4.248    0.0000   273.4   95.8
## 37 37 2005    10   2   12   122.8  318.5     0.000    2.6712   274.3   96.8
## 38 38 2005    10   2   13   104.4  325.9     0.000    2.5956   274.8   95.8
## 39 39 2005    10   2   14    63.3  322.8     0.000    3.1932   274.4   95.8
## 40 40 2005    10   2   15    46.4  324.3     0.000    2.4084   273.8   95.7
## 41 41 2005    10   2   16    25.8  317.9     0.000    3.2004   274.1   95.7
## 42 42 2005    10   2   17    12.2  323.0     0.000    2.8368   274.7   95.8
## 43 43 2005    10   2   18     0.3  326.4     0.000    1.1592   275.8   95.8
## 44 44 2005    10   2   19     0.0  328.6     0.000    1.5408   275.7   95.8
## 45 45 2005    10   2   20     0.0  329.7     0.000    2.5308   275.3   95.8
## 46 46 2005    10   2   21     0.0  328.6     0.000    0.6588   275.6   95.8
## 47 47 2005    10   2   22     0.0  330.0     0.000    0.0000   275.6   95.8
## 48 48 2005    10   2   23     0.0  330.6     0.000    0.2268   275.6   95.8
##    WindSpeed SurfPress       SWE SnowAlbedo     ShortOut  LongOut  SensFlux
## 25       1.3     87070        NA         NA           NA       NA        NA
## 26       1.3     87050        NA         NA           NA       NA        NA
## 27       1.0     87010        NA         NA           NA       NA        NA
## 28       2.1     86980        NA         NA           NA       NA        NA
## 29       1.5     86950        NA         NA           NA       NA        NA
## 30       1.8     86960        NA         NA           NA       NA        NA
## 31       1.3     86950        NA         NA           NA       NA        NA
## 32       1.2     86970        NA         NA           NA       NA        NA
## 33       1.4     86980        NA         NA           NA       NA        NA
## 34       1.5     86970        NA         NA           NA       NA        NA
## 35       1.9     86990        NA         NA           NA       NA        NA
## 36       1.9     86970 4.2457159  0.9500000  -41.1350000 -315.637  1.725718
## 37       1.9     86920 4.3515649  0.9477556 -116.3843896 -315.637  7.938301
## 38       2.7     86910 3.8994986  0.9455224  -98.7125412 -315.637 13.828172
## 39       2.3     86890 3.5260823  0.9433004  -59.7109136 -315.637  9.552238
## 40       1.9     86890 3.2737728  0.9410894  -43.6665483 -315.637  4.486866
## 41       1.7     86890 3.0643143  0.9388895  -24.2233481 -315.637  6.206740
## 42       2.5     86890 2.6880580  0.9367005  -11.4277460 -315.637 12.417438
## 43       4.1     86890 1.9865034  0.9345224   -0.2803567 -315.637 29.062365
## 44       3.1     86910 1.3159252  0.9323552    0.0000000 -315.637 23.255057
## 45       3.5     86910 0.6610287  0.9301989    0.0000000 -315.637 21.195883
## 46       4.1     86890        NA         NA           NA       NA        NA
## 47       4.6     86900        NA         NA           NA       NA        NA
## 48       3.4     86910        NA         NA           NA       NA        NA
##       LatFlux ColdContent       Melt   Runoff   SnowDepth LiqWatCont   IceCont
## 25         NA          NA         NA       NA          NA         NA        NA
## 26         NA          NA         NA       NA          NA         NA        NA
## 27         NA          NA         NA       NA          NA         NA        NA
## 28         NA          NA         NA       NA          NA         NA        NA
## 29         NA          NA         NA       NA          NA         NA        NA
## 30         NA          NA         NA       NA          NA         NA        NA
## 31         NA          NA         NA       NA          NA         NA        NA
## 32         NA          NA         NA       NA          NA         NA        NA
## 33         NA          NA         NA       NA          NA         NA        NA
## 34         NA          NA         NA       NA          NA         NA        NA
## 35         NA          NA         NA       NA          NA         NA        NA
## 36 -1.7987422           0 0.03744966 0.000000 0.041648549 0.03744966 4.2082662
## 37  3.8052949           0 0.28658825 2.570183 0.041871428 0.42505480 3.9265101
## 38  7.0398523           0 0.47200139 3.056606 0.036821677 0.43605044 3.4634482
## 39  3.9392064           0 0.33281813 3.571618 0.032688916 0.39045008 3.1356322
## 40  0.2261266           0 0.21478352 2.660997 0.029809376 0.35263695 2.9211359
## 41  1.7381866           0 0.18662678 3.412066 0.027416347 0.32759800 2.7367163
## 42  6.0733456           0 0.36357316 3.220769 0.023640441 0.30720265 2.3808553
## 43 18.5137561           0 0.68901821 1.884264 0.017179446 0.27115675 1.7153466
## 44 14.5996749           0 0.61846496 2.229917 0.011194638 0.20050427 1.1154210
## 45 12.3960564           0 0.60329986 3.201438 0.005533609 0.13316663 0.5278621
## 46         NA          NA         NA       NA          NA         NA        NA
## 47         NA          NA         NA       NA          NA         NA        NA
## 48         NA          NA         NA       NA          NA         NA        NA
##    SnowTemp SnowDensity            DateTime Precip
## 25       NA          NA 2005-10-02 00:00:00 0.4032
## 26       NA          NA 2005-10-02 01:00:00 0.0000
## 27       NA          NA 2005-10-02 02:00:00 0.0000
## 28       NA          NA 2005-10-02 03:00:00 0.0000
## 29       NA          NA 2005-10-02 04:00:00 0.0000
## 30       NA          NA 2005-10-02 05:00:00 0.7020
## 31       NA          NA 2005-10-02 06:00:00 0.9936
## 32       NA          NA 2005-10-02 07:00:00 3.1392
## 33       NA          NA 2005-10-02 08:00:00 2.4084
## 34       NA          NA 2005-10-02 09:00:00 1.4688
## 35       NA          NA 2005-10-02 10:00:00 3.4128
## 36   273.15    101.9415 2005-10-02 11:00:00 4.2480
## 37   273.15    103.9268 2005-10-02 12:00:00 2.6712
## 38   273.15    105.9023 2005-10-02 13:00:00 2.5956
## 39   273.15    107.8678 2005-10-02 14:00:00 3.1932
## 40   273.15    109.8236 2005-10-02 15:00:00 2.4084
## 41   273.15    111.7696 2005-10-02 16:00:00 3.2004
## 42   273.15    113.7059 2005-10-02 17:00:00 2.8368
## 43   273.15    115.6326 2005-10-02 18:00:00 1.1592
## 44   273.15    117.5496 2005-10-02 19:00:00 1.5408
## 45   273.15    119.4571 2005-10-02 20:00:00 2.5308
## 46       NA          NA 2005-10-02 21:00:00 0.6588
## 47       NA          NA 2005-10-02 22:00:00 0.0000
## 48       NA          NA 2005-10-02 23:00:00 0.2268

Das hat gut geklappt, aber was bedeutet die fehlende Angabe nach dem Komma?

Unser DataFrame ist in Zeilen und Spalten organisiert, die erste Dimension beinhaltet die Zeilen, hier haben wir mit 25:48 einen Bereich gewählt, die zweite Dimension die Spalten, hier haben wir keine Angabe gemacht, somit wurden alle Spalten gewählt!

Im Beispiel haben wir alle Stunden des 02.10.2005 ausgeschnitten! Aber natürlich muss ich hierzu vorher den zugehörigen Indexbereich kennen.

Um unabhängig von den Indizes ein Zeitfenster auszuwählen kann ich auch den Zeitraum definieren. Das Datum setzen uns wieder mit der paste-Funktion aus year, month und day zusammen und transformieren das Ergebnis dann gleich wieder in ein Datum:

# Datumsspalte erzeugen
input_data$Date = paste(input_data$Year,'-',input_data$Month,'-',input_data$Day,sep='')
input_data$Date = as.Date(input_data$Date)

# Bereich über Datum ausschneiden und Header ausgeben
time_subset = input_data[input_data$Date>"2005-10-01" & input_data$Date<"2005-10-03",]
head(time_subset, n=50)
##     X Year Month Day Hour ShortIn LongIn SolPrecip LiqPrecip AirTemp RelHum
## 25 25 2005    10   2    0     0.0  339.7     0.000    0.4032   277.5   97.1
## 26 26 2005    10   2    1     0.0  322.8     0.000    0.0000   277.1   97.0
## 27 27 2005    10   2    2     0.0  331.4     0.000    0.0000   277.1   97.0
## 28 28 2005    10   2    3     0.0  327.2     0.000    0.0000   276.9   94.9
## 29 29 2005    10   2    4     0.0  326.7     0.000    0.0000   276.7   93.9
## 30 30 2005    10   2    5     0.0  328.6     0.000    0.7020   276.1   88.8
## 31 31 2005    10   2    6     0.3  325.0     0.000    0.9936   276.1   91.8
## 32 32 2005    10   2    7     3.9  328.2     0.000    3.1392   274.7   94.8
## 33 33 2005    10   2    8    14.4  317.7     0.000    2.4084   273.8   96.8
## 34 34 2005    10   2    9    35.3  324.4     0.000    1.4688   274.9   97.9
## 35 35 2005    10   2   10    53.1  327.2     0.000    3.4128   274.2   96.9
## 36 36 2005    10   2   11    43.3  314.4     4.248    0.0000   273.4   95.8
## 37 37 2005    10   2   12   122.8  318.5     0.000    2.6712   274.3   96.8
## 38 38 2005    10   2   13   104.4  325.9     0.000    2.5956   274.8   95.8
## 39 39 2005    10   2   14    63.3  322.8     0.000    3.1932   274.4   95.8
## 40 40 2005    10   2   15    46.4  324.3     0.000    2.4084   273.8   95.7
## 41 41 2005    10   2   16    25.8  317.9     0.000    3.2004   274.1   95.7
## 42 42 2005    10   2   17    12.2  323.0     0.000    2.8368   274.7   95.8
## 43 43 2005    10   2   18     0.3  326.4     0.000    1.1592   275.8   95.8
## 44 44 2005    10   2   19     0.0  328.6     0.000    1.5408   275.7   95.8
## 45 45 2005    10   2   20     0.0  329.7     0.000    2.5308   275.3   95.8
## 46 46 2005    10   2   21     0.0  328.6     0.000    0.6588   275.6   95.8
## 47 47 2005    10   2   22     0.0  330.0     0.000    0.0000   275.6   95.8
## 48 48 2005    10   2   23     0.0  330.6     0.000    0.2268   275.6   95.8
##    WindSpeed SurfPress       SWE SnowAlbedo     ShortOut  LongOut  SensFlux
## 25       1.3     87070        NA         NA           NA       NA        NA
## 26       1.3     87050        NA         NA           NA       NA        NA
## 27       1.0     87010        NA         NA           NA       NA        NA
## 28       2.1     86980        NA         NA           NA       NA        NA
## 29       1.5     86950        NA         NA           NA       NA        NA
## 30       1.8     86960        NA         NA           NA       NA        NA
## 31       1.3     86950        NA         NA           NA       NA        NA
## 32       1.2     86970        NA         NA           NA       NA        NA
## 33       1.4     86980        NA         NA           NA       NA        NA
## 34       1.5     86970        NA         NA           NA       NA        NA
## 35       1.9     86990        NA         NA           NA       NA        NA
## 36       1.9     86970 4.2457159  0.9500000  -41.1350000 -315.637  1.725718
## 37       1.9     86920 4.3515649  0.9477556 -116.3843896 -315.637  7.938301
## 38       2.7     86910 3.8994986  0.9455224  -98.7125412 -315.637 13.828172
## 39       2.3     86890 3.5260823  0.9433004  -59.7109136 -315.637  9.552238
## 40       1.9     86890 3.2737728  0.9410894  -43.6665483 -315.637  4.486866
## 41       1.7     86890 3.0643143  0.9388895  -24.2233481 -315.637  6.206740
## 42       2.5     86890 2.6880580  0.9367005  -11.4277460 -315.637 12.417438
## 43       4.1     86890 1.9865034  0.9345224   -0.2803567 -315.637 29.062365
## 44       3.1     86910 1.3159252  0.9323552    0.0000000 -315.637 23.255057
## 45       3.5     86910 0.6610287  0.9301989    0.0000000 -315.637 21.195883
## 46       4.1     86890        NA         NA           NA       NA        NA
## 47       4.6     86900        NA         NA           NA       NA        NA
## 48       3.4     86910        NA         NA           NA       NA        NA
##       LatFlux ColdContent       Melt   Runoff   SnowDepth LiqWatCont   IceCont
## 25         NA          NA         NA       NA          NA         NA        NA
## 26         NA          NA         NA       NA          NA         NA        NA
## 27         NA          NA         NA       NA          NA         NA        NA
## 28         NA          NA         NA       NA          NA         NA        NA
## 29         NA          NA         NA       NA          NA         NA        NA
## 30         NA          NA         NA       NA          NA         NA        NA
## 31         NA          NA         NA       NA          NA         NA        NA
## 32         NA          NA         NA       NA          NA         NA        NA
## 33         NA          NA         NA       NA          NA         NA        NA
## 34         NA          NA         NA       NA          NA         NA        NA
## 35         NA          NA         NA       NA          NA         NA        NA
## 36 -1.7987422           0 0.03744966 0.000000 0.041648549 0.03744966 4.2082662
## 37  3.8052949           0 0.28658825 2.570183 0.041871428 0.42505480 3.9265101
## 38  7.0398523           0 0.47200139 3.056606 0.036821677 0.43605044 3.4634482
## 39  3.9392064           0 0.33281813 3.571618 0.032688916 0.39045008 3.1356322
## 40  0.2261266           0 0.21478352 2.660997 0.029809376 0.35263695 2.9211359
## 41  1.7381866           0 0.18662678 3.412066 0.027416347 0.32759800 2.7367163
## 42  6.0733456           0 0.36357316 3.220769 0.023640441 0.30720265 2.3808553
## 43 18.5137561           0 0.68901821 1.884264 0.017179446 0.27115675 1.7153466
## 44 14.5996749           0 0.61846496 2.229917 0.011194638 0.20050427 1.1154210
## 45 12.3960564           0 0.60329986 3.201438 0.005533609 0.13316663 0.5278621
## 46         NA          NA         NA       NA          NA         NA        NA
## 47         NA          NA         NA       NA          NA         NA        NA
## 48         NA          NA         NA       NA          NA         NA        NA
##    SnowTemp SnowDensity            DateTime Precip       Date
## 25       NA          NA 2005-10-02 00:00:00 0.4032 2005-10-02
## 26       NA          NA 2005-10-02 01:00:00 0.0000 2005-10-02
## 27       NA          NA 2005-10-02 02:00:00 0.0000 2005-10-02
## 28       NA          NA 2005-10-02 03:00:00 0.0000 2005-10-02
## 29       NA          NA 2005-10-02 04:00:00 0.0000 2005-10-02
## 30       NA          NA 2005-10-02 05:00:00 0.7020 2005-10-02
## 31       NA          NA 2005-10-02 06:00:00 0.9936 2005-10-02
## 32       NA          NA 2005-10-02 07:00:00 3.1392 2005-10-02
## 33       NA          NA 2005-10-02 08:00:00 2.4084 2005-10-02
## 34       NA          NA 2005-10-02 09:00:00 1.4688 2005-10-02
## 35       NA          NA 2005-10-02 10:00:00 3.4128 2005-10-02
## 36   273.15    101.9415 2005-10-02 11:00:00 4.2480 2005-10-02
## 37   273.15    103.9268 2005-10-02 12:00:00 2.6712 2005-10-02
## 38   273.15    105.9023 2005-10-02 13:00:00 2.5956 2005-10-02
## 39   273.15    107.8678 2005-10-02 14:00:00 3.1932 2005-10-02
## 40   273.15    109.8236 2005-10-02 15:00:00 2.4084 2005-10-02
## 41   273.15    111.7696 2005-10-02 16:00:00 3.2004 2005-10-02
## 42   273.15    113.7059 2005-10-02 17:00:00 2.8368 2005-10-02
## 43   273.15    115.6326 2005-10-02 18:00:00 1.1592 2005-10-02
## 44   273.15    117.5496 2005-10-02 19:00:00 1.5408 2005-10-02
## 45   273.15    119.4571 2005-10-02 20:00:00 2.5308 2005-10-02
## 46       NA          NA 2005-10-02 21:00:00 0.6588 2005-10-02
## 47       NA          NA 2005-10-02 22:00:00 0.0000 2005-10-02
## 48       NA          NA 2005-10-02 23:00:00 0.2268 2005-10-02

Um nun bestimmte Variablen aus dem DataFrame auszuwählen, kann ich diese bei der Auswahl wie folgt ansprechen (hier am Beispiel des Schneewasseräquivalents SWE), unsere Zeitinformation DateTime nehmen wir gleich mit:

# Bereich über Spaltennamen auswählen
variable_subset = input_data[,c("DateTime","SWE")]
head(variable_subset, n=50)
##               DateTime       SWE
## 1  2005-10-01 00:00:00        NA
## 2  2005-10-01 01:00:00        NA
## 3  2005-10-01 02:00:00        NA
## 4  2005-10-01 03:00:00        NA
## 5  2005-10-01 04:00:00        NA
## 6  2005-10-01 05:00:00        NA
## 7  2005-10-01 06:00:00        NA
## 8  2005-10-01 07:00:00        NA
## 9  2005-10-01 08:00:00        NA
## 10 2005-10-01 09:00:00        NA
## 11 2005-10-01 10:00:00        NA
## 12 2005-10-01 11:00:00        NA
## 13 2005-10-01 12:00:00        NA
## 14 2005-10-01 13:00:00        NA
## 15 2005-10-01 14:00:00        NA
## 16 2005-10-01 15:00:00        NA
## 17 2005-10-01 16:00:00        NA
## 18 2005-10-01 17:00:00        NA
## 19 2005-10-01 18:00:00        NA
## 20 2005-10-01 19:00:00        NA
## 21 2005-10-01 20:00:00        NA
## 22 2005-10-01 21:00:00        NA
## 23 2005-10-01 22:00:00        NA
## 24 2005-10-01 23:00:00        NA
## 25 2005-10-02 00:00:00        NA
## 26 2005-10-02 01:00:00        NA
## 27 2005-10-02 02:00:00        NA
## 28 2005-10-02 03:00:00        NA
## 29 2005-10-02 04:00:00        NA
## 30 2005-10-02 05:00:00        NA
## 31 2005-10-02 06:00:00        NA
## 32 2005-10-02 07:00:00        NA
## 33 2005-10-02 08:00:00        NA
## 34 2005-10-02 09:00:00        NA
## 35 2005-10-02 10:00:00        NA
## 36 2005-10-02 11:00:00 4.2457159
## 37 2005-10-02 12:00:00 4.3515649
## 38 2005-10-02 13:00:00 3.8994986
## 39 2005-10-02 14:00:00 3.5260823
## 40 2005-10-02 15:00:00 3.2737728
## 41 2005-10-02 16:00:00 3.0643143
## 42 2005-10-02 17:00:00 2.6880580
## 43 2005-10-02 18:00:00 1.9865034
## 44 2005-10-02 19:00:00 1.3159252
## 45 2005-10-02 20:00:00 0.6610287
## 46 2005-10-02 21:00:00        NA
## 47 2005-10-02 22:00:00        NA
## 48 2005-10-02 23:00:00        NA
## 49 2005-10-03 00:00:00        NA
## 50 2005-10-03 01:00:00        NA

6 - Daten aggregieren

Oft liegen Daten zeitlich feiner aufgelöst (z.B. als Stundenwerte) vor als wir sie bräuchten (z.B. Monatswerte). Um die Daten zeitlich zu aggregieren, braucht es ein Zeitkriterium anhand dessen aggregiert werden kann, z.B. eine Spalte die nur Jahre- und Monatsinformation enthält - ich kann dann z.B. den Mittelwert oder die Summe über dieses Kriterium berechnen.

Wir testen das anhand unserer Niederschlagsdaten, dafür fügen wir unserem DataFrame eine Spalte hinzu, die nur Jahr und Monat enthält und verwenden die Funktion yearmon() des Paketes zoo um diese Zeichenfolge in eine Zeitinformation zu transformieren:

# Bibliothek einbinden
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
# Spalte mit Jahr und Monat erstellen
input_data$YearMonth = paste(input_data$Year,'-',input_data$Month,sep='')
input_data$YearMonth = as.yearmon(input_data$YearMonth,format="%Y-%m")

# Header ausgeben
head(input_data$YearMonth, n=10)
##  [1] "Oct 2005" "Oct 2005" "Oct 2005" "Oct 2005" "Oct 2005" "Oct 2005"
##  [7] "Oct 2005" "Oct 2005" "Oct 2005" "Oct 2005"

Nun verwenden wir die aggregate() Funktion um die Niederschlagssumme für die Monate zu berechnen, das Ergebnis schauen wir gleich mit head() an:

# Monatliche Niederschläge berechnen
monthly_Precip <- aggregate(input_data$Precip,FUN = sum,by = list(Group.date = input_data$YearMonth))

# Header ausgeben
head(monthly_Precip, n=10)
##   Group.date         x
## 1   Oct 2005 164.82636
## 2   Nov 2005  99.27468
## 3   Dec 2005 158.77692
## 4   Jan 2006  97.14528
## 5   Feb 2006 109.91376
## 6   Mar 2006 136.73599
## 7   Apr 2006  26.43555
## 8   May 2006  49.80496
## 9   Jun 2006  52.51841

Die Spaltenüberschriften passen wir noch so an, dass man auch weiss was hier enthalten ist! Wir verwenden hierzu die Funktion setNames():

# Spaltennamen definieren
monthly_Precip = setNames(monthly_Precip,c("Month","Precip"))

# Header ausgeben
head(monthly_Precip, n=10)
##      Month    Precip
## 1 Oct 2005 164.82636
## 2 Nov 2005  99.27468
## 3 Dec 2005 158.77692
## 4 Jan 2006  97.14528
## 5 Feb 2006 109.91376
## 6 Mar 2006 136.73599
## 7 Apr 2006  26.43555
## 8 May 2006  49.80496
## 9 Jun 2006  52.51841

Die hier erhaltenen Monatswerte stellen wir nun als Balkendiagramm in blau dar. Verwendet dabei den Typ type="h" sowie das Argument lwd um die Dicke Eurer Balken anzupassen:

Tipp:

Die Farbe kann wie gewohnt einfach unter Angabe einer Farbe gesetz werden (col="blue") oder über einen RGB Wert definiert werden:

col=rgb(0,0,0,alpha=0.3)

Diese Option bietet auch die Möglichkeit einen Transparenzwert alpha zu definieren, allerdings ist aufzupassen, da RGB Werte in der Regel einen Bereich von 0-255 haben, hier aber eine Intensität von 0-1 eingesetzt werden muss! Man muss die RGB-Werte falls im 0-255er Wertebereich vorliegend, also normieren, d.h. durch den Maximalwert teilen:

Für die Farbe blau heisst das: rgb(0,0,255) –> rgb(0/255,0/255,255/255)

Eine Übersicht über verschiedene Farben und deren RGB Werte findet Ihr unter http://www.farb-tabelle.de/de/farbtabelle.htm!

# Monatlichen Niederschlag plotten
plot(monthly_Precip$Month, monthly_Precip$Precip, type="h", col=rgb(0/255,0/255,255/255, alpha=0.5), ann=FALSE, lwd=20)

# Diagrammtitel hinzufügen
title(main="Monatlicher Niederschlag \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Monat")
title(ylab="Niederschlag (mm)")

Der Plot sieht soweit ganz gut aus, lediglich die Zeitachse ist noch eine Mischung aus Monat und Jahr - R tut das um uns den Jahreswechsel zu zeigen, aber das sieht nicht so schön aus. Wir werden die Formatierung der Zeitachse also selbst in die Hand nehmen, dafür geben wir in unserem plot-Aufruf den Parameter xaxt="n" mit so dass zunächst keine Achsenbeschriftung erfolgt, dann fügen wir am Schluss ein eigenes Axen-Element mit entsprechend formatierter Zeitinformation über den Befehl axis() hinzu:

# Monatlichen Niederschlag plotten
plot(monthly_Precip$Month, monthly_Precip$Precip, type="h", col=rgb(0/255,0/255,255/255, alpha=0.5), xaxt = "n", ann=FALSE, lwd=20)

# Diagrammtitel hinzufügen
title(main="Monatlicher Niederschlag \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Monat")
title(ylab="Niederschlag (mm)")

# Add dates to x-axis
axis(1,monthly_Precip$Month,format(monthly_Precip$Month, "%Y-%m"))

7 - Scatterplots

Scatterplots (deutsch: Streudiagramme) ermöglichen es den Zusammenhang zweier Größen darzustellen. Dabei wird die abhängige Variable (Y) als Funktion der unabhängigen Variable (X) dargestellt.

Wir wollen nun Scatterplots auf Basis simulierter und beobachteter Schneegrößen erstellen - diese Art der Auswertung stellt eine wichtige Säule in der Evaluation von Modellergebnissen dar.

Wir lesen dazu beobachte Werte des Schneewasseräquivalents zum Vergleich mit den Simulationen ein - diese Daten befinden sich in der Datei Observations.txt im Arbeitsverzeichnis.

Schaut Euch die Daten zunächst im Editor an und versucht den Lesebefehl richtig zu konfigurieren, ohne unten nachzusehen!

# Beobachtungsdaten einlesen
observation_file = "Observations.txt"
daily_observed_SWE = read.table(observation_file, sep="", dec = ".", header=TRUE)

Wir sehen uns das Ergebnis des Einlesens mit head() an:

# Header ausgeben
head(daily_observed_SWE, n=10)
##    year month day SWE
## 1  2005    10   1   0
## 2  2005    10   2   0
## 3  2005    10   3   0
## 4  2005    10   4   0
## 5  2005    10   5   0
## 6  2005    10   6   0
## 7  2005    10   7   0
## 8  2005    10   8   0
## 9  2005    10   9   0
## 10 2005    10  10   0

Wir sehen dass alles geklappt hat, allerdings liegen die Beobachtungen als Tageswerte vor - wir müssen also zunächts die Simulationen auf Tageswerte aggregieren. Dafür verwenden wir wieder die aggregate() Funktion zusammen mit unserem vorhin erzeugten Datum:

# Stündliche SWE Simulationen auf Tageswerte aggregieren 
daily_simulated_SWE <- aggregate(input_data$SWE,FUN = mean,by = list(Group.date = input_data$Date))

Ein Blick auf das Ergebnis:

# Header ausgeben
head(daily_simulated_SWE, n=60)
##    Group.date        x
## 1  2005-10-01       NA
## 2  2005-10-02       NA
## 3  2005-10-03       NA
## 4  2005-10-04       NA
## 5  2005-10-05       NA
## 6  2005-10-06       NA
## 7  2005-10-07       NA
## 8  2005-10-08       NA
## 9  2005-10-09       NA
## 10 2005-10-10       NA
## 11 2005-10-11       NA
## 12 2005-10-12       NA
## 13 2005-10-13       NA
## 14 2005-10-14       NA
## 15 2005-10-15       NA
## 16 2005-10-16       NA
## 17 2005-10-17       NA
## 18 2005-10-18       NA
## 19 2005-10-19       NA
## 20 2005-10-20       NA
## 21 2005-10-21       NA
## 22 2005-10-22       NA
## 23 2005-10-23       NA
## 24 2005-10-24       NA
## 25 2005-10-25       NA
## 26 2005-10-26       NA
## 27 2005-10-27       NA
## 28 2005-10-28       NA
## 29 2005-10-29       NA
## 30 2005-10-30       NA
## 31 2005-10-31       NA
## 32 2005-11-01       NA
## 33 2005-11-02       NA
## 34 2005-11-03       NA
## 35 2005-11-04       NA
## 36 2005-11-05       NA
## 37 2005-11-06       NA
## 38 2005-11-07       NA
## 39 2005-11-08       NA
## 40 2005-11-09       NA
## 41 2005-11-10       NA
## 42 2005-11-11       NA
## 43 2005-11-12       NA
## 44 2005-11-13       NA
## 45 2005-11-14       NA
## 46 2005-11-15       NA
## 47 2005-11-16       NA
## 48 2005-11-17       NA
## 49 2005-11-18       NA
## 50 2005-11-19       NA
## 51 2005-11-20       NA
## 52 2005-11-21       NA
## 53 2005-11-22       NA
## 54 2005-11-23       NA
## 55 2005-11-24       NA
## 56 2005-11-25 20.19373
## 57 2005-11-26 25.36564
## 58 2005-11-27 26.39569
## 59 2005-11-28 28.41766
## 60 2005-11-29 33.74083

Das Ergebnis schaut gut aus, wieder sind die Spaltenbeschriftungen nicht sehr aussagekräftig, wir ändern diese deshalb unter Verwendung der setNames() Funktion:

# Spaltenbeschriftung anpassen
daily_simulated_SWE = setNames(daily_simulated_SWE,c("Date","SWE"))

Wir schauen kurz auf die Anzahl der Messungen und Simulationen, in unserem Fall sollten diese beiden Zeitreihen gleich lag sein, wäre das nicht so müssten wir einen Überlapppungszeitraum verwenden! Wir verwenden hierzu die Funktion nrow(), welche die Anzahl der Zeilen eines DataFrames ausgibt:

# Anzahl der Zeilen ausgeben
nrow(daily_simulated_SWE)
## [1] 273
nrow(daily_observed_SWE)
## [1] 273

Die beiden Zeitreihen sind gleich lang (auch wenn das nicht unbedingt heisst, dass die Tage die gleichen sind), wir erstellen uns nun eine neue Liste mit den SWE Beobachtungen und Simulationen, nur der Übersichtlichkeit wegen:

# Liste mit Beobachtungen und Simulationen anlegen
scatter_input = NULL
scatter_input$Observed_SWE = daily_observed_SWE$SWE
scatter_input$Simulated_SWE = daily_simulated_SWE$SWE

# Struktur ausgeben
str(scatter_input)
## List of 2
##  $ Observed_SWE : num [1:273] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Simulated_SWE: num [1:273] NA NA NA NA NA NA NA NA NA NA ...

Diese Liste übergeben wir nun an die plot() Funktion, dabei plotten wir die Beobachtungen auf die X-Achse und die Simulationen auf die Y-Achse. Dabei können wir Vieles anwenden, das wir eben über das Formatieren von Plots gelernt haben:

# Scatterplot erzeugen
plot(scatter_input$Observed_SWE, scatter_input$Simulated_SWE, type="p", col="black", ann=FALSE)

# Diagrammtitel hinzufügen
title(main="Simuliertes vs. beobachtetes Schneewasseräquivalent \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Beobachtetes SWE (mm)")
title(ylab="Simuliertes SWE (mm)")

# Grid hinzufügen
grid(lty = 'dotted', col = "gray", lwd = 1)

Unser Scatterplot schaut schon ganz gut aus, ein paar Details wollen wir aber noch ändern:

  1. Bei einem Scatterplot muss der Wertebereich der X- und Y-Achse immer gleich sein, nur so wäre das perfekte Modell eine Linie im 45 Grad Winkel!

  2. Wir wollen den Wertebereich auf minimal 0 umstellen, negative SWE-Werte gibt es nicht und somit brauchen wir auch keinen Bereich < 0

  3. Die Kreissymbole wollen wir größer machen und auch füllen

Den Wertebereich können wir mit den bereits bekannten Parameter xlim bzw. ylim anpassen, damit ist sicher gestellt, dass X- und Y-Achse gleich skaliert sind.

Merke:

Wenn wir später unseren Plot speichern, stellen wir durch Wahl einer Höhe und Breite sicher, dass unser Plot bei gleicher X- und Y-Skalierung dann auch quadratisch ist, also dass das Verhältnis 1:1 sichergestellt ist. Wenn wir in RStudio selbst plotten, ist das Grafikfenster “Plots” je nach Rechner/Bildschirm unterschiedlich groß (wir können dieses ja sogar größer und kleiner ziehen, was den Plot neu skaliert). Wenn wir auch hier in der Polt-Ansicht in jedem Fall wollen, dass unser Plot in einem Achsenverhältnis 1:1 dargestellt wird, können wir das durch das Setzen von par(pty="s") for dem Plot Befehl tun. Alternativ gibt es einen Parameter asp=, wenn wir diesen im Plotaufruf auf 1 setzen, wird das Verhältnis zwischen X- und Y-Achse auf 1 gesetzt, aber dieser Parameter ignoriert unsere xlim bzw. ylim Einstellungen und ist somit hier nicht hilfreich.

Die Art des Symbols können wir über den Parameter pch einstellen. Folgende Optionen sind hier möglich (wollt Ihr diese Symbole als Grafik sehen oder mehr über den Umgang mit diesen erfahren besucht: https://r-charts.com/base-r/pch-symbols/):

Die Größe des Symbols kann über cex eingestellt werden, wobei der Wert 1 den Standard repräsentiert, größere (z.B. 1.5) oder kleinere (z.B. 0.5) Werte vergrößern/verkleinern entsprechend.

# Scatterplot erzeugen
par(pty="s")
plot(scatter_input$Observed_SWE, scatter_input$Simulated_SWE, xlim=c(0,600), ylim=c(0,600), type="p", pch = 19, cex=1.2, col="black",ann=FALSE)

# Diagrammtitel hinzufügen
title(main="Simuliertes vs. beobachtetes Schneewasseräquivalent \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Beobachtetes SWE (mm)")
title(ylab="Simuliertes SWE (mm)")

# Grid hinzufügen
grid(nx = NULL, ny = NULL, lty = 'dotted', col = "gray", lwd = 1)

Um den Zusammenhang zwischen Simulationen und Beobachtungen statistisch zu beschreiben, braucht es (hier) eine lineare Regressionsfunktion. Diese wollen wir uns im Folgenden näher ansehen.

8 - Lineare Regression

Die lineare Regression versucht einen Y-Wert in linearer Abhängigkeit eines X-Wertes zu berechnen.

Heraus kommt eine Regressionsfunktion:

y = Achsenabschnitt + Steigung * x = b + m * x

Solche Regressionsfunktionen kommen zum Beispiel bei Streudiagrammen (wie hier) oder bei Trendanalysen in Klimawandelstudien zum Einsatz. In Letzteren versucht man damit die Veränderung einer abhängigen Variablen (z.B. der Temperatur, Y-Achse) in Abhängigkeit einer unabhängigen Variablen (den Jahren, X-Achse) zu untersuchen.

Um die lineare Beziehung zwischen Simulationen und Beobachtungen herzuleiten verwenden wir die Funktion lm() (linear model) und übergeben dieser zunächst die Y-Werte (die Simulationen) und dann die X-Werte (die Beobachtungen). Die Ergebnisse schreiben wir in die Variable fit_model:

# Lineare Regression berechnen
fit_model = lm(scatter_input$Simulated_SWE ~ scatter_input$Observed_SWE)

Durch diesen Aufruf haben wir eine Variable fit_model definiert, die alle Information zur berechneten Regression enthält. Die Koeffizienten (Achsenabschnitt = Intercept) und die Steigung (hier: scatter_input$Observed_SWE) kann man sich über Aufruf von coefficients() und Angabe der Regressionsfunktion fit_model ausgeben lassen:

# Koeffizieneten ausgeben
coefficients(fit_model)
##                (Intercept) scatter_input$Observed_SWE 
##                  20.759101                   1.121946

Was heisst das nun? Das bedeutet, ausgehend von einem Wert von 20 (mm SWE) an der Y-Achse steigt unsere Regressionsgerade mit einer Steigung von 1.12 mit jedem weiteren beobachteten mm SWE! Aber schauen wir uns das später gleich im Plot genauer an…

Will man sich gleich die gesamte Zusammenfassung der Linearen Regression ansehen, verwendet man die Funktion summary() unter Angabe der Regressionsfunktion fit_model.

Ausgeben werden dann:

  1. unter Call die verwendete Funktion
  2. die Statistik der Abweichungen (Residuen) unter Residuals
  3. die Regressionskoeffizienten (Achsenabschnitt und Steigung) zusammen mit den Ergebnissen des Signifikanztests unter Coefficients
  4. der Standardfehler, das R2 und die F-Statistik.

Aber was heisst Signifikanz eigentlich: Der hier verwendete Signifikanztest testet die Regressionskoeffizienten gegen “0”. Er testet also, mit welcher Wahrscheinlichkeit ich richtig liege dass z.B. die Steigung ungleich “0” ist, also ein deutlicher Zusammenhang besteht.

Wir schauen uns nun die Zusammenfassung der Regression einmal an:

# Zusammenfassung ausgeben
summary(fit_model) 
## 
## Call:
## lm(formula = scatter_input$Simulated_SWE ~ scatter_input$Observed_SWE)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.062 -27.555  -9.564   8.808 106.406 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                20.75910    7.18493   2.889  0.00441 ** 
## scatter_input$Observed_SWE  1.12195    0.02784  40.299  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39.73 on 157 degrees of freedom
##   (114 observations deleted due to missingness)
## Multiple R-squared:  0.9118, Adjusted R-squared:  0.9113 
## F-statistic:  1624 on 1 and 157 DF,  p-value: < 2.2e-16

Wir sehen, dass ein signifikanter Zusammenhang besteht, das Signifikanzniveau ist mit mind. zwei Sternen angegeben, also min. 0.001. Wir gehen also mit einer Irrtumswahrscheinlichkeit von 0.1 % fälschlicherweise von einem Zusammenhang der beiden Größen aus, wir liegen also mit einer Wahrscheinlichkeit von 99.9% richtig wenn wir einen Zusammenhang annehmen!

Wir wollen nun die Regressionslinie auch im Plot einfügen und ansehen. Dazu plotten nochmal unser Streudiagramm und verwenden die Funktion abline() um die Regressionslinie hinzuzufügen. Wir geben dieser einen gestrichelten Linientyp lty='dashed', eine Linienstärke von 3 lwd=3 und die Farbe rot col='red'.

# Scatterplot erzeugen
par(pty="s")
plot(scatter_input$Observed_SWE, scatter_input$Simulated_SWE, xlim=c(0,600), ylim=c(0,600), type="p", pch = 19, cex=1.2, col="black", ann=FALSE)

# Diagrammtitel hinzufügen
title(main="Simuliertes vs. beobachtetes Schneewasseräquivalent \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Beobachtetes SWE (mm)")
title(ylab="Simuliertes SWE (mm)")

# Grid hinzufügen
grid(lty = 'dotted', col = "gray", lwd = 1)

# Regressionslinie hinzufügen
abline(fit_model, lty = 'dashed', col = 'red', lwd = 3)

Zu einer Regressionsgerade gehört aber immer die zugehörige Regressionsfunktion, und am besten auch das Bestimmtheitsmaß (R2), beides wollen wir nun als Text hinzufügen.

Dazu müssen wir zunächst durch Zuhilfenahme der Funktion round(Variable,Dezimalstellen) die in der Variable coefficients gespeicherten Regressionskoeffizienten sowie das unter r.squared zu findende Bestimmtheitsmaß auf zwei Stellen hinter dem Komma runden.

Danach basteln wir uns eine Textvariable regress_function in die wir mit der Funktion paste() nacheinander Textbausteine und Variablen einbauen.

Die Funktion text() erlaubt dann das Einfügen der generierte Zeichenfolge in das Diagramm.

Nun das Ganze Schritt für Schritt:

Zunächst runden wir alle Variablen, stellen den Text zusammen und geben diesen Testweise ohne Diagramm aus:

# Parameter der Regressionsgleichung auf 2 Dezimalstellen runden
slope = round(summary(fit_model)$coefficients[2], 2)
intercept = round(summary(fit_model)$coefficients[1], 2)
rsquared = round(summary(fit_model)$r.squared, 2)

# Zeichenfolge mit Regressionsgleichung definieren
regress_function = paste("y = ", intercept, " + ", slope, " * x",sep='')

# Zeichenfolge mit Bestimmtheitsmaß definieren
rsquared_text = paste("R2 = ", rsquared, sep='')

# Text ausgeben
print(regress_function)
## [1] "y = 20.76 + 1.12 * x"
print(rsquared_text)
## [1] "R2 = 0.91"

Nun bauen wir diesen Text in unser Diagramm ein. Dies tun wir unter Verwendung der text() Funktion, diese benötigt die X-Koordinaten im Plot, die Y-Koordinaten im Plot und den Text in Form von:

text(X-Koordimnate,Y-Koordinate,MeinText)

Durch Hinzufügen des Parameters pos können wir noch die Position definieren, d.h. ob der Text unter 1, links 2, über 3 oder rechts 4 von den angegeben Koordinaten platziert werden soll. Wir vewenden den Wert 4 (rechts).

Wollt Ihr mehr über die text-Funktion erfahren, besucht: https://www.rdocumentation.org/packages/graphics/versions/3.6.2/topics/text

Die X- und Y-Koordinaten sind dabei über dei Werte der Achsen zu skalieren, sie stellen also Werte in der Plotfläche dar!

Nach Abgleich mit unseren Achsenwerten setzen wir die Regressionsfunktion auf die X-koordinate 0 und die Y-Koordinate 550. Den Text mit dem R2 setzen wir knapp darunter auf die X-koordinate 0 und die Y-Koordinate 500. Bis auf die text-Funktion ganz unten verwenden wir unseren vorherigen plot-Aufruf!

Tipp:

Während die Angabe der Koordinaten nach Sichtung des Achsenbereichs gut über absolute Werte klappt, geht das natürlich nicht mehr wenn sich der Wertebereich ändert! Man kann mit etwas Aufwand die Koordinaten relativ zum Wertebereich setzen, natürlich nur wenn man diesen vorher über xlim bzw. ylim definiert, also kennt. Wenn ich meinen Text z.B. immer auf 95% des Wertebereiches (der Y-Achse) setzen will kann ich das so tun:

ypos = ymin + (0.95 * (ymax - ymin))

Analog funktiomiert das auch mit der X-Achse!

# Scatterplot erzeugen
par(pty="s")
plot(scatter_input$Observed_SWE, scatter_input$Simulated_SWE, xlim=c(0,600), ylim=c(0,600), type="p", pch = 19, cex=1.2, col="black", ann=FALSE)

# Diagrammtitel hinzufügen
title(main="Simuliertes vs. beobachtetes Schneewasseräquivalent \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Beobachtetes SWE (mm)")
title(ylab="Simuliertes SWE (mm)")

# Grid hinzufügen
grid(nx = NULL, ny = NULL, lty = 'dotted', col = "gray", lwd = 1)

# Regressionslinie hinzufügen
abline(fit_model, lty = 'dashed', col = 'red', lwd = 3)

# Regressionsfunktion als Text hinzufügen
text(0, 550, regress_function, pos=4)

# Bestimmtheitsmaß als Text hinzufügen
text(0, 520, rsquared_text, pos=4)

Das hat gut geklappt! Allerdings ist das Bestimmtheitsmaß allein nur wenig geeignet um Modellergebnisse zu bewerten, wir wollen im Folgenden deshalb weitere Effizienzkriterien kennenlernen!

9 - Effizienzkriterien

Effizienzkriterien können herangezogen werden, um Modellergebnisse zu bewerten. Folgende Kriterien sollten wir kennen:

  1. Bestimmtheitsmaß R2, Wertebereich: 0 <–> 1:

    R2 (das Quadrat des Korrelationskoeffizienten R) sagt uns welcher Teil der Gesamtvarianz durch die abhängige Variable (im Plot: Y-Achse) erklärt wird. In unserem Fall erklärt das Modell 91% der Gesamtvarianz. R2 hat aber große Nachteile, da es nur die Kovarianz betrachtet, aber nie die absoluten Differenzen zwischen Beobachtung und Modell. So kann das R2 1 sein, wenn ein Modell immer 1000 mm SWE über den Beobachtungen liegt, da beide auch dann perfekt kovarieren, mit einer guten Modellperformanz hat das natürlich dann nichts zu tun.

  2. Nash-Suttcliffe Koeffizient NSE, Wertebereich: minus undendlich <–> 1:

    Der Nash-Sutcliffe Koeffizient (auch Nash-Sutcliffe Model Efficiency, NSE) betrachtet auch die Differenzen zwischen Beobachtung und Modellergebnis und ist somit kritischer und in der Hydrologie ein Standard Kriterium. Das Perfekte Modell hätte einen NSE Wert von 1, ab einem Wert von 0 ist der Mittelwert eine bessere Vorhersage als das Modell.

  3. PBIAS, Wertebereich: minus unendlich <–> plus unendlich:

    Der PBIAS gibt den Grad der Über-/Unterschätzung durch das Modell an. Ein PBIAS von -20 beschreibt eine Unterschätzung von 20%, ein PBIAS von +20 beschreibt eine Überschätzung von 20%. Das perfekte Modell liefert einen PBIAS Wert von 0.

Das waren nur einige wichtige Vertreter, für einen umfassenden Überblick zu verschiedenen Effizienzkriterien in der Hydrologie, schaut Euch folgende Artikel an:

https://adgeo.copernicus.org/articles/5/89/2005/

https://swat.tamu.edu/media/1312/moriasimodeleval.pdf

Um NSE zu berechnen könnt Ihr die Bibliothek topmodel verwenden. Die Formel zur Berechnung von NSE lautet:

\[ NSE=1 - [\sum_{i = 1}^{n}{(obs_i-sim_i)^2}: \sum_{i = 1}^{n}{(obs_i-obs_m)^2}] \] Dabei stehen obs und sim für die Beobachtungen und Simulationen, i und m für die Werte eines Zeitschrittes bzw. das Mittel aller Werte.

Verwendet doch gleich einmal die enthaltene Funktion NSeff(Qobs,Qsim) um NSE zu berechnen und setzt NSE noch unter dem R2 in Euren Plot ein!

Zur Berechnung des PBIAS habe ich Euch eine kleine Funktion gebaut, die mit PBIAS(Qobs,Qsim) gerufen werden kann. Sie tut nichts anderes, als für jeden Zeitschritt die Abweichung des Simulationswertes vom Beobachtungswert zu berechnen und diese Abweichungen dann aufzusummieren. Die Summe der Abweichungen wird dann geteilt durch die Summe aller Messwerte, mit 100 Multipliziert ergibt sich die durchschnittliche prozentuale Abweichung:

\[ PBIAS=100 * [\sum_{i = 1}^{n}{(sim_i-obs_i)} : \sum_{i = 1}^{n}{(obs_i)}] \]

Die unten aufgeführte Funktion müsst Ihr in Eurem Skript einbauen, nachdem man die Definition einmal ausgeführt hat, kann man die Funktion immer wieder aufrufen.

# define function
PBIAS = function(Qobs,Qsim) {

# get indices with valid values
valid_obs=which(!is.na(Qobs))
valid_sim=which(!is.na(Qsim))
valid_values=intersect(valid_obs,valid_sim)

# get bias
bias=Qsim[valid_values]-Qobs[valid_values]
sum_bias=sum(bias)
sum_obs=sum(Qobs[valid_values])
pbias=100*sum_bias/sum_obs

# return result
return(pbias)
}

# use funtion to derive PBIAS
pbias=PBIAS(scatter_input$Observed_SWE,scatter_input$Simulated_SWE)
pbias
## [1] 21.14462

Wir überschätzen das beobachtete SWE also im Schnitt um 21%. Verwendet nun diese Funktion um Euren PBIAS zu berechnen und setzt den berechneten Wert noch unter dem NSE in Euren Plot ein!

Eine grobe Richtlinie zur Bewertung von Modellergebnissen über NSE und PBIAS geben Morasi et al. (2007):

10 - Plots speichern

Für das Speichern von Plots gibt es verschiedene Möglichkeiten.

  1. Wir können das händisch in RStudio tun, indem wir in der Plotansicht auf Export->Save as Image gehen

  2. Wenn wir im Jupyter arbeiten, könnt Ihr die generierten Grafiken auch mit Rechtsklick --> Speichern unter ... im Browser speichern

  3. Alternativ (und für Skripte meist besser weil automatisert möglich) können wir den Plot im Code selbst durch den png() Befehl (vor dem Plotten einzufügen) gefolgt von dev.off() (nach dem Plotten einzufügen) speichern

Option 3) probieren wir gleich mal anhand eines Plots aus, wir verwenden dafür nochmal den oben generierten Scatterplot, geben dabei der png() Funktion auch noch die Breite und Höhe unseres pngs mit, da die Qualität (Auflösung) standardmäßig redzuiert ist. Beachtet dass wir durch eine geeignete Wahl der Breite und Höhe dafür sorgen müssen, dass unser Achsenverhältnis nicht verzerrrt ist, wir wählen eine leicht erhöhte Höhe um dem Titel Platz zu geben:

# Datei aufmachen
png(file = "Scatterplot.png", width=500, height=520)

# Scatterplot erzeugen
par(pty="s")
plot(scatter_input$Observed_SWE, scatter_input$Simulated_SWE, xlim=c(0,600), ylim=c(0,600), type="p", pch = 19, cex=1.2, col="black", ann=FALSE)

# Diagrammtitel hinzufügen
title(main="Simuliertes vs. beobachtetes Schneewasseräquivalent \nStation Col de Porte 2005/2006")

# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Beobachtetes SWE (mm)")
title(ylab="Simuliertes SWE (mm)")

# Grid hinzufügen
grid(nx = NULL, ny = NULL, lty = 'dotted', col = "gray", lwd = 1)

# Regressionslinie hinzufügen
abline(fit_model, lty = 'dashed', col = 'red', lwd = 3)

# Regressionsfunktion als Text hinzufügen
text(0, 550, regress_function, pos=4)

# Bestimmtheitsmaß als Text hinzufügen
text(0, 520, rsquared_text, pos=4)

# NSE als Text hinzufügen
library(topmodel)
nse_text=paste("NSE = ", round(NSeff(scatter_input$Observed_SWE,scatter_input$Simulated_SWE),2), sep='')
text(0, 490, nse_text, pos=4)

# PBIAS als Text hinzufügen
pbias_text=paste("PBIAS = ", round(PBIAS(scatter_input$Observed_SWE,scatter_input$Simulated_SWE),2), sep='')
text(0, 460, pbias_text, pos=4)

# Device schliessen
dev.off()
## quartz_off_screen 
##                 2

Die so erzeugte Datei Scatterplot.png liegt nun im Arbeitsverzeichnis und sollte etwa so aussehen: